
library(tidyverse)
library(modelsummary)
library(broom)
library(marginaleffects)
library(psych)
library(ggeffects)
library(devtools)
library(estimatr)


#Insert data without respondents who failed the attention check.
#Those who dplyr::selected an option other than 5 for question Q20.1.6.

data <- read.csv("replication_ssq_isekihandhata.csv") %>% 
  filter(Q20.1_6 == 5)


table(data$Q20.1_6)


#Variables related to the level of trust in experts/governing party/opposition party.
data <- data %>% 
  mutate(Experts = if_else(Q9.1_1==6 |Q9.1_1 ==7, NA_real_, 6 - Q9.1_1)) %>% 
  mutate(Trust_LDP = if_else(Q9.1_2 ==6 |Q9.1_2 ==7, NA_real_, 6 - Q9.1_2)) %>% 
  mutate(Trust_CDP = if_else(Q9.1_3 ==6 |Q9.1_3 ==7, NA_real_, 6 - Q9.1_3))


#Variables related to the necessity of the state of emergency declaration.
data <- data %>% 
  mutate(unnecessary = if_else(Q4.2_6==6 |Q4.2_6 ==7, NA_real_, 6 - Q4.2_6))

#Covariate
data <- data %>% 
  mutate(edu = if_else(Q2.2 == 6, NA_real_, Q2.2),
         income = as.factor(Q2.3),
         cs = if_else(Q2.4 == 5, NA_real_, Q2.4))

#Party Identification
data <- data %>% 
  mutate(pid = 
           case_when(Q7.1 == 1 | Q7.1 == 5 ~ "The governing \n party supporter",
                     Q7.1 == 10 | Q7.1 == 11 ~ NA_character_,
                     Q7.1 ==  9 ~ "Independent",
                     TRUE ~ "The opposition \n party supporter"))

###########################################################################################################################
###########################################################################################################################
#Regression Analysis on Trust in Experts and Political Parties

fit1 <- data %>% 
  lm_robust(Experts ~ pid * unnecessary + sex + age  + cs + edu + income, 
            data = ., 
            se_type = "stata")

fit2 <- data %>% 
  lm_robust(Trust_LDP ~ pid * unnecessary + sex + age  + cs + edu + income, 
            data = ., 
            se_type = "stata")

fit3 <- data %>% 
  lm_robust(Trust_CDP ~ pid * unnecessary + sex + age  + cs + edu + income, 
            data = ., 
            se_type = "stata")

#Estimate Marginal Effect
me1 <- fit1 %>% 
  plot_slopes(., "unnecessary", condition = "pid", draw = TRUE) 

me2 <- fit2 %>% 
  plot_slopes(., "unnecessary", condition = "pid", draw = TRUE) 

me3 <- fit3 %>% 
  plot_slopes(., "unnecessary", condition = "pid", draw = TRUE) 


result.Experts <- tibble(me1$plot_env$dat) %>% 
  mutate(team = "Trust in Experts") %>% 
  dplyr::select(estimate,conf.low,conf.high, pid, team)

result.LDP <- tibble(me2$plot_env$dat) %>% 
  mutate(team = "Trust in LDP") %>% 
  dplyr::select(estimate,conf.low,conf.high, pid, team)

result.CDP <- tibble(me3$plot_env$dat) %>% 
  mutate(team = "Trust in CDP") %>% 
  dplyr::select(estimate,conf.low,conf.high, pid, team)


#Combine the data frames of each fitting value
result_ALL <- rbind(result.Experts,result.LDP,result.CDP) %>% 
  mutate(sig = if_else(conf.low * conf.high > 0, "sig", "ns")) %>% 
  mutate(team = factor(team, levels = c("Trust in Experts",
                                        "Trust in LDP",
                                        "Trust in CDP"))) %>% 
  mutate(pid = factor(pid, levels = c("The governing \n party supporter",
                                      "The opposition \n party supporter",
                                      "Independent")))

#Draw a figure
figure1 <- result_ALL %>% 
  ggplot() +
  geom_hline(yintercept = 0, color = "black",linetype="dashed") + 
  geom_pointrange(aes(x = pid, ymin = conf.low, ymax = conf.high,
                      y = estimate, color = sig),) +
  scale_color_manual(values = c("sig" = "red",
                                "ns" = "black")) +
  theme_gray(base_size = 16) +
  theme(legend.position = "bottom") + 
  facet_wrap( ~ team) + 
  labs(x = "Party Identification",
       y = "Marginal Effect") + 
  theme_bw() + 
  guides(color = F)

figure1

ggsave("Figure1.png", figure1, dpi = 300, width = 13, height = 5)


#############################################################
#Estimate the predicted value using {marginaleffects}
#note: The figures generated by the code below differ slightly from those 
# appearing in the published paper due to a version upgrade of the {marginaleffects} package.

pre1 <- predictions(fit1,
                    newdata = datagrid(pid = 
                                         c("The governing \n party supporter",
                                           "The opposition \n party supporter",
                                           "Independent"), 
                                       unnecessary = c(1:5))) %>% 
  mutate(outcome = "Trust in Experts") %>% 
  dplyr::select(estimate,unnecessary,conf.low,conf.high,pid,outcome)

pre2 <- predictions(fit2,
                    newdata = datagrid(pid = 
                                         c("The governing \n party supporter",
                                           "The opposition \n party supporter",
                                           "Independent"), 
                                       unnecessary = c(1:5))) %>% 
  mutate(outcome = "Trust in LDP") %>% 
  dplyr::select(estimate,unnecessary,conf.low,conf.high,pid,outcome)

pre3 <- predictions(fit3,
                    newdata = datagrid(pid = 
                                         c("The governing \n party supporter",
                                           "The opposition \n party supporter",
                                           "Independent"), 
                                       unnecessary = c(1:5))) %>% 
  mutate(outcome = "Trust in CDP") %>% 
  dplyr::select(estimate,unnecessary,conf.low,conf.high,pid,outcome)


#Combine the data frames of each predicted value
pre <- rbind(pre1,pre2,pre3) %>% 
  mutate(outcome = factor(outcome, levels = c("Trust in Experts",
                                        "Trust in LDP",
                                        "Trust in CDP"))) %>% 
  mutate(pid = factor(pid, levels = c("The governing \n party supporter",
                                      "The opposition \n party supporter",
                                      "Independent")))

#Draw a figure
figure2 <- pre %>% 
  ggplot() +
  geom_line(aes(y = estimate, x = unnecessary)) + 
  geom_ribbon(aes(y = estimate, ymin = conf.low, ymax = conf.high,
                  x = unnecessary),alpha=0.3) +
  geom_text(aes(x = unnecessary, y = estimate, 
                label = sprintf("%2.2f", estimate)), size = 3.5,                                                          
            position = position_stack(vjust = 0)) +
  # scale_y_continuous(breaks = seq(1,4,4),limits=c(1,4)) +
  labs(x = "I don't think declaring an emergency is necessary.", y = "Predicted Value") + 
  facet_grid(pid ~ outcome) + 
  theme_bw(base_size = 12)

figure2
ggsave("Figure2.png", figure2, dpi = 300, width = 10, height = 7)



###########################################################################################################################
###########################################################################################################################
#
